Get started

## Open libraries
## general basics
library(here)
library(ggplot2)

## the basic packages for spatial analyses and viz
library(sf)
library(terra)
library(tmap)

## for telemetry and home range analyses
library(move)
library(adehabitatHR) 
library(sp)

Define workspace and folder structure:

work_wd <- here()
data_raw_wd <- here("data-raw")
geo_raw_wd  <- here("data-raw", "geo-raw")
output_wd   <- here("output")

Load data sets:

# if you are not sure about the file names, check the folder content with
# list.files(data_raw_wd)
animals <- st_read(layer = "DieTiere_32633", dsn = data_raw_wd)  

## Reading layer `DieTiere_32633' from data source 
##   `C:\Users\DataVizard\Google Drive\Work\Eco\Projects\2022_BioDivCourse\Course4_MoveQ\data-raw' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 388 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 390397.5 ymin: 5862042 xmax: 396748.2 ymax: 5868366
## Projected CRS: WGS 84 / UTM zone 33N

landcov <- st_read(layer = "Landnutzung_3035", dsn = geo_raw_wd)    

## Reading layer `Landnutzung_3035' from data source 
##   `C:\Users\DataVizard\Google Drive\Work\Eco\Projects\2022_BioDivCourse\Course4_MoveQ\data-raw\geo-raw' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 182 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 4539362 ymin: 3304861 xmax: 4566087 ymax: 3324860
## Projected CRS: ETRS89-extended / LAEA Europe

rivers  <- st_read(layer = "linearegew_32633", dsn = geo_raw_wd)   

## Reading layer `linearegew_32633' from data source 
##   `C:\Users\DataVizard\Google Drive\Work\Eco\Projects\2022_BioDivCourse\Course4_MoveQ\data-raw\geo-raw' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2834 features and 13 fields
## Geometry type: MULTILINESTRING
## Dimension:     XYZ
## Bounding box:  xmin: 381989.2 ymin: 5852103 xmax: 408856.6 ymax: 5871804
## z_range:       zmin: 0 zmax: 0
## Projected CRS: WGS 84 / UTM zone 33N

Define or transform coordinate reference system (CRS):

lc_32633 <- st_transform(landcov, 32633)

Make a plot of the data:

## Create a color code for the 5 land use classes
land_cols <- c("red", "antiquewhite", "chartreuse4", "cornsilk4", "blue")

## TODO - argh, Punkte und Linien erscheinen nicht im html.... 
## CED: Hab leider keinen Schimmer von base pltos... ggplot for the rescue?

plot(lc_32633[,"NS1"], col = land_cols[lc_32633$NS1], border = "transparent")
plot(animals[,1], col = animals$ID, add = TRUE) 
plot(rivers[,1], col = 'blue', lwd = 0.1, add = TRUE)

ggplot2:

## Create a color code for the 5 land use classes
land_cols <- c("red", "antiquewhite", "chartreuse4", "cornsilk4", "navyblue")

## Create a color code for the 3 individuals
id_cols <- viridis::inferno(n = 12)[c(2,4,6)]

ggplot() +
  geom_sf(data = lc_32633, aes(fill = as.factor(NS1)), color = NA, alpha = .4) +
  geom_sf(data = rivers, color = "navyblue", alpha = .4) +
  geom_sf(data = animals, aes(color = as.factor(ID))) +
  scale_color_manual(values = id_cols, name = "Animal ID") +
  scale_fill_manual(values = land_cols, name = "Land Use Class") +
  theme_minimal(base_size = 15)

tmap:

## Create a color code for the 5 land use classes
land_cols <- c("red", "antiquewhite", "chartreuse4", "cornsilk4", "navyblue")

## Create a color code for the 3 individuals
id_cols <- viridis::inferno(n = 12)[c(2,4,6)]

tmap_mode(mode = "view")

tm_shape(shp = lc_32633[,'NS1']) + 
    tm_polygons(col = "NS1", palette = land_cols, alpha = .4) +
  tm_shape(shp = rivers[,1]) + 
    tm_lines(col = 'blue')  +
  tm_shape(shp = animals[,'ID']) + 
    tm_dots(col = "ID", palette = id_cols)

• a. Calculate their home ranges.
  
  and
  
  
• b. Buffer the home range with the mean displacement per step of all three animals.
  

First , we retrieve the step length. For this, transform the sf object into a move object:

anim_coords <- as.data.frame(st_coordinates(animals)) ## returns a list
animals$X <- anim_coords$X
animals$Y <- anim_coords$Y
animals_ng <- st_drop_geometry(animals)

animals_move <- move(x = animals_ng$X, y = animals_ng$Y,
                     time = as.POSIXct(animals_ng$MYDATUM, 
                                       format = "%Y-%m-%d 8:0:0", tz = "CET"), ## CED: Puh, das ist aber komplex... geht's nicht iwie einfacher, z.B. mit lubridate?
                     proj = CRS("+init=epsg:32633"),
                     data = animals_ng, animal = animals_ng$ID, sensor = "GPS")

Check the data by plotting the new move object:

plot(animals_move, type = "l")

Since we have 3 animals, a list with 3 elements is returned:

steplength_per_animal <- distance(animals_move)
mean_stepl_p_anim <- lapply(steplength_per_animal,mean)
(mean(unlist(mean_stepl_p_anim)))

## [1] 303.0478

The mean step length is approx. 300 m.

Now use this value to buffer the MCPs:

## for using mcp() in {adehabitatHR}, convert to SpatialPointsDataframe of {sp}
animals_sp <- as(animals, "Spatial")

## calculate MCP
mcp <- mcp(animals_sp[,"ID"], percent = 95) 

## the buffer command needs however an sf object
mcp_sf  <- as(mcp, "sf")
mcp_buf <- st_buffer(mcp_sf, 300)

## the plot is correct
plot(mcp)                         
plot(animals_sp, add = TRUE, col = animals_sp$ID) 
plot(mcp_buf[,1], add = TRUE)
plot(animals_sp, add = TRUE, col = id_cols) 
plot(mcp, add = TRUE)

## save buffer to disc
st_write(obj = mcp_buf, 
         dsn = output_wd, 
         layer = "mcp_buffer_32633",
         driver = "ESRI Shapefile", 
         delete_layer = TRUE)

## Deleting layer `mcp_buffer_32633' using driver `ESRI Shapefile'
## Writing layer `mcp_buffer_32633' to data source 
##   `C:/Users/DataVizard/Google Drive/Work/Eco/Projects/2022_BioDivCourse/Course4_MoveQ/output' using driver `ESRI Shapefile'
## Writing 3 features with 2 fields and geometry type Polygon.

• c. Clip the home ranges with the underlying landscape and calculate the percentage of land use classes 2 (agriculture) and 2 (forest) within the home range (hint: use st_intersection()).
  

## check the total area of the buffer in km2 per animal:
buf_area <- st_area(mcp_buf) / 1e6

## clip with land use
mcp_buf_lc <- st_intersection(lc_32633, mcp_buf) 

## have a look at the object:
mcp_buf_lc

## Simple feature collection with 45 features and 14 fields
## Geometry type: GEOMETRY
## Dimension:     XY
## Bounding box:  xmin: 390158.2 ymin: 5861742 xmax: 397048.2 ymax: 5868666
## Projected CRS: WGS 84 / UTM zone 33N
## First 10 features:
##        AREA PERIMETER CLCGES_BES CLCGES_B_1  NS BRD SZENE   SPHEROID   TKNR NS1
## 2  33397200  79888.10      23163      23162 312   1 19323 Krassovsky N33111   3
## 5  11562800  44001.70      23265      23264 211   1 19323 Krassovsky N33111   2
## 15 18208500  73976.60      23764      23763 231   1 19323 Krassovsky N33111   2
## 16  8568800  37160.80      23788      23787 231   1 19323 Krassovsky N33111   2
## 28  5142930  13701.10      24028      24027 312   1 19323 Krassovsky N33111   3
## 35  1227710   4956.00      24314      24313 243   1 19323 Krassovsky N33111   2
## 36  1448920  11624.70      24348      24347 311   1 19323 Krassovsky N33111   3
## 41   654073   3560.28      24431      24430 211   1 19323 Krassovsky N33111   2
## 42  3534710  11425.00      24454      24453 311   1 19323 Krassovsky N33111   3
## 44 46474200 125933.00      24480      24479 211   1 19323 Krassovsky N33111   2
##    NS2 NS3 id     area                       geometry
## 2   31 312  1 1444.323 MULTIPOLYGON (((394196.9 58...
## 5   21 211  1 1444.323 POLYGON ((392577.5 5866547,...
## 15  23 231  1 1444.323 POLYGON ((391413.6 5866716,...
## 16  23 231  1 1444.323 POLYGON ((395139.3 5866196,...
## 28  31 312  1 1444.323 POLYGON ((390833 5866283, 3...
## 35  24 243  1 1444.323 POLYGON ((392292.4 5866469,...
## 36  31 311  1 1444.323 POLYGON ((395809.1 5865790,...
## 41  21 211  1 1444.323 POLYGON ((394729.2 5865833,...
## 42  31 311  1 1444.323 POLYGON ((390395.2 5865319,...
## 44  21 211  1 1444.323 MULTIPOLYGON (((392247.2 58...

## take care, you have to recalculate the area!
mcp_buf_lc$areakm2 <- st_area(mcp_buf_lc) / 1e6

## now, for each animal id, summarize the area per NS1-category 2 and 3:
agri   <- tapply(mcp_buf_lc$areakm2[mcp_buf_lc$NS1 == 2], ## CED: Das können wir wohl von denen nicht erwarten oder??
                 mcp_buf_lc$id[mcp_buf_lc$NS1 == 2], 
                 FUN = sum)

forest <- tapply(mcp_buf_lc$areakm2[mcp_buf_lc$NS1 == 3],
                 mcp_buf_lc$id[mcp_buf_lc$NS1 == 3], 
                 FUN = sum)

## percentages
round(agri / buf_area * 100, 0)

## Units: [1/m^2]
##  1  2  3 
## 71 48 85

round(forest / buf_area * 100, 0)

## Units: [1/m^2]
##  1  2  3 
## 29 52 14

• d. In which land use are the animals actually located? Extract land use underneath each telemetry location (hint: st_join()). Do the animals actually avoid agriculture?
  

spatjoin <- st_join(animals, lc_32633)
myres <- table(spatjoin$NS1, spatjoin$ID)
myres

##    
##       1   2   3
##   2  48  25  95
##   3 147  64   9

## now divide that by the total number of locations per animal
myn <- table(animals$ID)

round((myres[1,] / myn*100), 0)  ## agriculture

## 
##  1  2  3 
## 25 28 91

round((myres[2,] / myn*100), 0)  ## forest

## 
##  1  2  3 
## 75 72  9

Compare the values with the percentages above.

For example, animal 1 had 71% of agriculture and 29% of forest in its HR, however, only 25% of the locations were actually in agriculture, but 75% of its locations were in forest. We could conclude that animal 1 avoids agriculture.

Sys.time()

## [1] "2022-03-07 18:21:20 CET"

sessionInfo()

## R version 4.1.2 (2021-11-01)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19043)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252   
## [3] LC_MONETARY=German_Germany.1252 LC_NUMERIC=C                   
## [5] LC_TIME=C                      
## system code page: 65001
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] adehabitatHR_0.4.19 adehabitatLT_0.3.25 CircStats_0.2-6    
##  [4] boot_1.3-28         MASS_7.3-54         adehabitatMA_0.3.14
##  [7] ade4_1.7-18         deldir_1.0-6        move_4.1.6         
## [10] rgdal_1.5-28        raster_3.5-11       sp_1.4-6           
## [13] geosphere_1.5-14    tmap_3.3-2          terra_1.4-22       
## [16] sf_1.0-5            ggplot2_3.3.5       here_1.0.1         
## 
## loaded via a namespace (and not attached):
##  [1] RColorBrewer_1.1-2      httr_1.4.2              rprojroot_2.0.2        
##  [4] tools_4.1.2             bslib_0.3.1             utf8_1.2.2             
##  [7] R6_2.5.1                KernSmooth_2.23-20      DBI_1.1.2              
## [10] colorspace_2.0-2        withr_2.4.3             gridExtra_2.3          
## [13] tidyselect_1.1.1        leaflet_2.0.4.1         compiler_4.1.2         
## [16] leafem_0.1.6            textshaping_0.3.6       xml2_1.3.2             
## [19] bookdown_0.24           sass_0.4.0              scales_1.1.1           
## [22] classInt_0.4-3          proxy_0.4-26            systemfonts_1.0.3      
## [25] stringr_1.4.0           digest_0.6.29           rmarkdown_2.11         
## [28] base64enc_0.1-3         dichromat_2.0-0         pkgconfig_2.0.3        
## [31] htmltools_0.5.2         highr_0.9               fastmap_1.1.0          
## [34] htmlwidgets_1.5.4       rlang_0.4.12            farver_2.1.0           
## [37] jquerylib_0.1.4         generics_0.1.1          jsonlite_1.7.2         
## [40] crosstalk_1.2.0         dplyr_1.0.7             magrittr_2.0.1         
## [43] s2_1.0.7                Rcpp_1.0.7              munsell_0.5.0          
## [46] fansi_0.5.0             viridis_0.6.2           abind_1.4-5            
## [49] lifecycle_1.0.1         stringi_1.7.5           leafsync_0.1.0         
## [52] yaml_2.2.1              tmaptools_3.1-1         grid_4.1.2             
## [55] parallel_4.1.2          crayon_1.4.2            lattice_0.20-45        
## [58] stars_0.5-5             knitr_1.36              pillar_1.6.4           
## [61] codetools_0.2-18        wk_0.6.0                XML_3.99-0.8           
## [64] glue_1.4.2              evaluate_0.14           leaflet.providers_1.9.0
## [67] png_0.1-7               vctrs_0.3.8             rmdformats_1.0.3       
## [70] gtable_0.3.0            purrr_0.3.4             assertthat_0.2.1       
## [73] cachem_1.0.6            xfun_0.27               lwgeom_0.2-8           
## [76] e1071_1.7-9             ragg_1.1.3              class_7.3-19           
## [79] viridisLite_0.4.0       tibble_3.1.6            memoise_2.0.1          
## [82] units_0.7-2             ellipsis_0.3.2